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Abstract 

Background: Type 1 diabetes (T1D) is an autoimmune disease, while type 2 (T2D) and gestational diabetes (GDM) 
are considered metabolic disturbances. In a previous study evaluating the transcript profiling of peripheral 
mononuclear blood cells obtained from T1D, T2D and GDM patients we showed that the gene profile of T1D 
patients was closer to GDM than to T2D. To understand the influence of demographical, clinical, laboratory, 
pathogenetic and treatment features on the diabetes transcript profiling, we performed an analysis integrating 
these features with the gene expression profiles of the annotated genes included in databases containing 
information regarding GWAS and immune cell expression signatures. 

Methods: Samples from 56 (19 T1D, 20 T2D, and 17 GDM) patients were hybridized to whole genome one-color 
Agilent 4x44k microarrays. Non-informative genes were filtered by partitioning, and differentially expressed genes 
were obtained by rank product analysis. Functional analyses were carried out using the DAVID database, and 
module maps were constructed using the Genomica tool. 

Results: The functional analyses were able to discriminate between T1D and GDM patients based on genes 
involved in inflammation. Module maps of differentially expressed genes revealed that modulated genes: i) 
exhibited transcription profiles typical of macrophage and dendritic cells; ii) had been previously associated with 
diabetic complications by association and by meta-analysis studies, and iii) were influenced by disease duration, 
obesity, number of gestations, glucose serum levels and the use of medications, such as metformin. 

Conclusion: This is the first module map study to show the influence of epidemiological, clinical, laboratory, 
immunopathogenic and treatment features on the transcription profiles of T1D, T2D and GDM patients. 

Keywords: Type 1 diabetes, Type 2 diabetes, Gestational diabetes, Module map, Immune cell signatures, 
Transcriptome analysis 
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Background 

Diabetes mellitus is considered by the World Health 
Organization to be a global epidemic. Although the different 
types of diabetes have varying incidence rates, it is estimated 
that approximately 6.4% of the world adult population have 
some type of diabetes [1]. Diabetes is characterized as hyper- 
glycemia resulting from a relative or absolute impairment of 
insulin secretion as well as peripheral resistance to insulin 
[2]. Based on etiology, diabetes mellitus has been classified 
into type 1 diabetes (T1D), type 2 diabetes (T2D), gesta- 
tional diabetes mellitus (GDM), as well as other types of 
diabetes, including genetic defects in (3-cell function, gen- 
etic defects in insulin action, diseases of the exocrine pan- 
creas, endocrinopathies, drug- or chemical-induced forms, 
infection-induced diabetes, uncommon forms of immune- 
mediated diabetes, and other genetic syndromes associ- 
ated with diabetes [2]. The pathogenic mechanisms of 
each diabetes type are still unclear, especially in the case 
of T2D and GDM. Pregnant women with gestational dia- 
betes also have an increased risk of developing T2D, sug- 
gesting a close relationship between these types [3,4]. 

T1D is the most studied type of diabetes, has several 
susceptibility loci identified in mice (Iddl-Idd26) and is 
characterized by the autoimmune destruction of pancre- 
atic beta cells leading to insulin deficiency [5]. Macro- 
phages, dendritic cells and lymphocytes are involved in 
this pathogenic process through a complex interplay of 
mechanisms implicated in the loss of immune tolerance 
to autoantigens, including i) the hypoexpression of insu- 
lin in the thymus during promiscuous antigen expres- 
sion [5,6]; ii) autoantigen presentation mediated by 
molecules coded by the HLA-DRB1 *04-DQB 1:03:02 and 
HLA-DR3r03-DQB 1*02:01 haplotypes resulting in the 
development of insulin autoantibodies (IAA) and auto- 
antibodies against the 65 kDa isoform of glutamic acid 
decarboxylase (GAD A), respectively [7]; iii) a deregulation 
of the immune response mediated by either an impaired 
expression of surface regulatory molecules (IL2RA, IL2RB, 
and CTLA-4) or a deregulation of intracellular signals 
(PTPN2 and PTPN22) [8]; iv) a decreased number of sup- 
pressive or T regulatory cells [9]; v) a decreased number of 
iNKT cells [10]; and vi) a loss of function of molecules in- 
volved in the innate immune response [11,12], In addition, 
several other genes have been implicated in the develop- 
ment of T1D by human genome wide association studies 
(GWAS) [13-15], including genes identified by transcrip- 
tome analyses in human or animal models [16,17] that eval- 
uated peripheral blood mononuclear cells [18,19], 
pancreatic beta cells [20] and whole blood cells [21]. 

In recent years, genetic mechanisms have also been 
shown to affect the development of T2D [22]. Several 
studies have associated polymorphisms in the PPARG 
and KCNJ11 genes with T2D susceptibility [23], Other 
studies have identified a strong effect of TCF7L2 gene 



variants on T2D risk, possibly affecting proglucagon ex- 
pression with consequent reduced insulin secretion [23,24]. 
Indeed, both GWAS and international collaborative efforts 
to analyze GWAS data from multiple groups, such as the 
Meta- Analysis of Glucose and Insulin-related traits Consor- 
tium (MAGIC), have identified other genetic variants asso- 
ciated with T2D gene susceptibility [25,26], several of which 
were associated with glycemic traits. Many of these groups 
of genes were related to abnormal insulin processing 
(MADD), higher proinsulin and lower insulin secretion 
(TCF7L2, SLC30A8, GIPR, and C2CD4B), and abnor- 
malities in early insulin secretion (MTNR1B, FADS1, 
DGKB, and GCK) [22]. Although GWAS have identified 
susceptibility regions across the genome and transcriptome 
studies have indicated several modulated genes in beta cells 
and the usefulness of blood RNA profiles [27,28], functional 
studies are still needed to understand the role of several 
genes obtained from association studies. 

The pathogenesis of T1D has been considered to be 
different from that of T2D, and information obtained by 
GWAS has indicated that most T1D and T2D genetic 
loci seem to not overlap. However, there is evidence that 
inflammatory processes involving interleukin-1 may play 
a role in islet beta cell loss in both types of diabetes [29]. 
Peripheral blood mononuclear cell (PBMC) energy 
metabolism has been found to play a major role in the 
pathogenesis of insulin resistance [30]. In addition, 
there is increasing evidence that metabolic regulation in 
these cells could influence the number of PBMCs, prolif- 
eration pathways, molecular basal synthesis and leukocyte 
function [31]. In accordance with these findings, transcrip- 
tome studies comparing T1D with T2D [32] and T1D with 
a control group [18] have revealed changes in the expres- 
sion of several genes related to inflammatory response, 
fatty acid biosynthesis, hydrolase activity, detoxification of 
aldehydes generated by alcohol metabolism and lipid per- 
oxidation, all of which could affect the metabolism of 
PBMCs in diabetic patients. 

Compared with T1D and T2D, GDM has been subject 
to fewer linkage studies and transcriptome analyses 
[33]. Many genes, however, are known to be associated 
with both GDM and T2D, especially those related to 
obesity and oxidative stress [34,35]. The mechanisms 
linking excess adiposity to an elevated risk of GDM are 
not completely understood, but recent evidence points 
to the role of specific hormones and cytokines known as 
adipokines, which are secreted by the adipose tissue 
[36,37]. In addition, transcriptome signatures obtained 
from placenta [38] and whole blood cells [39] have iden- 
tified genes involved with lipid metabolism that are dif- 
ferentially expressed between T1D and GDM. Recently, 
in a meta-analysis of the transcription profiles of T1D, T2D 
and GDM patients, our group reported that gene expres- 
sion signatures of GDM patients were closer to those of 
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T1D patients than to T2D [40]. The analyses of gene ex- 
pression signatures, however, were impaired by the presence 
of multiple variables associated with each type of diabetes. 
To circumvent these problems, here we used several bio- 
informatics tools to analyze demographical, clinical, labora- 
tory, pathogenetic and treatment data against modulated 
genes that have been annotated by databases containing in- 
formation regarding both GWAS and gene expression sig- 
natures displayed by immune cells. This type of comparison 
has yielded multiple informative modules that have been 
used to perform comprehensive maps in cancer [41], and 
here we apply these tools to diabetes. 

Results 

Overall, a schematic heatmap with all demographic, clinical 
and laboratory patient features is shown in Figure 1. The 
global partitioning analysis of the three types of diabetes 
disclosed 8,469 transcripts considered as informative, which 
are available at www.rge.fmrp.usp.br/passos/DBF-MCL, the 
principal component analysis (PCA) of these genes are 
shown in Figure 2 and the summarized DAVID functional 
categories (Kegg pathways) of the main clusters are shown 
in Figure 3. 

The Venn diagrams yielded shared and specific genes 
after statistical analysis by rank products (T1D versus 
T2D, T2D versus GDM and T1D versus GDM) (Figure 4) 
as well as multiple significant summarized DAVID func- 
tional categories (Kegg pathways) (Figure 5). The mod- 
ule maps encompassing all analyses, i.e., genes obtained 
from both partitioning and rank products, were created 
with the set of genes obtained in each of the approaches 
described above (Figure 6). Finally, the confirmation by 
PCR analysis of important genes involved in diabetes is 
shown in Figure 7. 

Global analysis of the data set with density-based filter 
and Markov clustering (DBF-MCL) and principal component 
analysis (PCA) 

A list of 8,469 genes obtained from global partitioning 
analysis (DBF-MCL method applying a non-informative 
filter) yielded informative non-supervised clusters of 
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Figure 2 Principal component analysis (PCA) of the major types 
of diabetic patients, using the 8,469 informative genes obtained 
by the DBF-MCL algorithm. The separation of samples of each type 
of diabetes after filtering non-informative genes showed similarities 
among them, indicating that sample transcription profiles were not 
influenced by the batch effect. 



co-regulated genes. Figures showing all non-supervised 
clusters and the respective gene list are available at 
www.rge.fmrp.usp.br/passos/DBF-MCL. A PCA plot of 
filtered genes is shown (Figure 2). Despite we have more 
similarity between T1D and GDM, some patients (with 
the lowest levels of glycemic levels) were closer to the 
other types. Overall, many gene expression profiles are 
shared between T1D and GDM. A typical example of this 
analysis can be observed in cluster #1 (3,089 probes), in 
which all T2D patients shared a group of induced genes 
encoding zinc finger proteins, while these same genes 
were found repressed in T1D and GDM patients. Simi- 
larly, in cluster #12 (216 probes), all T2D patients shared a 
group of repressed genes associated with cytokine and 
chemokine activity, NOD -like and Toll-like receptor sig- 
naling pathways and MAPK signaling pathway (Figure 3). 

Supervised analysis of T1D, T2D and GDM patients 

To understand the influence of demographic, clinical, 
pathogenetic and laboratory features on the differential 
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Figure 1 Heatmap representative of type of diabetes, demographic, clinical, laboratory and treatment features of the patients. 

Qualitative variables were assigned by the absence or presence of the characteristic, and quantitative variables were assigned by values below or 
above the mean values. This information was used as array (experimental) set for the module map construction. 
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:T cell receptor signaling pathway 
:Neuroactive ligand-receptor interaction 
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Figure 3 Heatmap of the significant functional categories of the genes of the clusters obtained by non-informative filters, in the basis 
of a FDR < 0.1. The Kegg categories were obtained from DAVID knowledge base with an enrichment P value < 0.05 after Benjamini correction. 
The grey scale represents the logarithm of the enriched P value. 



gene expression profiles among diabetes patients, a sec- 
ond strategy was used. Module maps were created by 
comparing the three groups at the same time (from the 
partitioning analysis) as well as by the individual two by 
two group comparisons. These paired analysis module 
maps were constructed based on the results obtained 
from statistical analyses of rank products and are repre- 
sented in Venn diagrams (Figure 4), disclosing up and 
down-regulated genes shared by the same type of dia- 
betes in different analysis. The functional categories of 
these genes are shown in Figure 5. 

T1D patients presented several genes that were induced 
when compared to GDM patients, including genes in the 
MHC region (HLA-DQA1 and HLA-DQA2), TNF receptors 
(TNFRSF17, TNFAIP6), cytokines (ILIA, IL-IB, IL1F10, IL4, 
116, US, IL23A, IL27) and cytokine receptors (IL1R2, 
IL1R1N, IL18R1), chemokines (CXCL1, CXCL2, CCL20, 
CCL23, CCL3L3, CCL4) and chemokine receptors (CCR3, 
CRL2), lymphocyte receptors (KIR3DL2, KIR2DS4), and 
transcription factors (GLIS2, SOX8, GATA2, RUNXI, 
SOD2, FOXC1, FOXC2, FOXE1). Additionally, T1D pa- 
tients presented several genes that were differentially 
expressed compared to T2D patients, including immune 



response genes, as MHC region genes (HLA-DQA1 and 
DQA2), clusters of differentiation genes (CD8B, CD55, 
CD83), cytokines (ILIA, IL1BJL6, IL8, IL23A), chemokines 
(CXCL1, CCL3L3, CXCL2, CXCL3, CCL20, CCL23, CCL24, 
CCL3, CCL4), NOD-like receptor signaling pathway 
(CXCLI, IL6, CARD9, TNF, IL8, CXCL2, CASP8, NFKBIA, 
IL1B, MAPK8, TNFAIP3, NLRP3) and transcription fac- 
tors (GLIS2, SOX8, FOXC2, FOXC1, FOXE1, FOXK1). 

GDM patients presented several genes that were induced 
compared to T1D patients, including toll-like receptors 
(TLR6, TLR7), carbohydrate binding genes (LGALS3, 
LGALSI2, CLEC7A, CLECIB), plasma membrane recep- 
tors (ILSRA, CCR1), MHC genes (HLA-DRB3), solute car- 
rier family genes (SLC6A10P, SLC6A4, SLCIAS, SLC4A1, 
SLC8A1, SLC6A17, SLC16A3, SLC6A8, SLC14A1), cluster 
of differentiation genes (CD9, CD33, CD36), and chemo- 
kines (CXCLS, CCLIS, CXCLI2, CXCLI0). In contrast, 
GDM patients presented several genes that were induced 
compared to T2D, including immune response genes 
(CXCLI0, CXCR4, CD46, CCL3L3, IL1B, ILIA, NFKBIZ, 
IL27, C1QA, CD83, C1QB, HIF1A, DEFA4, DEFA3, 
KIR2DL4, CXCLI, CCL3, TNF, CCR1, CXCL3, CXCL2, 
CCL4, CCL23, CCL20, IL6, CR2, IL8, IL1RN STAT3, 
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Figure 4 Venn diagrams show the differentially expressed genes after paired analysis of the three types of DM. The genes were 
identified by Rank Product analysis with P value < 0.001 and a percentage of false prediction (pfp) < 0.05. The analysis referring to upregulated 
genes is shown in panel A and that of downregulated genes in panel B. 
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hsa04060 : Cytokine-cytokine receptor interaction 
hsa04621 :N0D-like receptor signaling pathway 
hsa04940:Type I diabetes mellitus 
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hsa04512 :ECM-receptor interaction 
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hsa04062 :Chemokine signaling pathway 
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Figure 5 Heatmap of the significative functional categories of the differentially expressed genes obtained by paired Rank Products 
analysis with P value < 0.001 and percentage of false prediction (pfp) < 0.05 (T1 D vs. GDM; T2D vs. GDM and T1 D vs. T2D). The Kegg 
categories were obtained from DAVID knowledge base with an enrichment pvalue< 0.05 after Benjamini correction. The grey scale represents 
the logarithm of the enriched lvalue. 
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Figure 6 Heatmaps of the modules identified by Genomica tool, which compares gene lists of immune cells and diabetic association 
genes with demographic, clinical, laboratory and treatment features of patients (P value < 0.05, corrected by the false discovery 
rate - FDR < 0.05. The four module maps presented list of genes (induced or repressed), identified by: A) Non-informative filters (DBF-MCL 
algorithm) in the basis of FDR < 10%; B) Rank products analysis of T1D vs. GDM; C) Rank products analysis of T2D vs. GDM; D) Rank Products analysis 
of T1D vs. T2D. Abbreviations: MF - macrophages ; Bla and B1 b - subsets of B lymphocytes; BFo - follicular B lymphocytes; BMz - marginal zone B 
lymphocytes; Treg - regulatory T lymphocytes; DC - dendritic cells ; CD4 and CD8 - subsets of T lymphocytes. 
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Figure 7 Confirmation of microarray findings by qRT-PCR of (A) IL1B, (B) RGS1, (C) EG/?2, (D) FOX03A, (E) SOD2 and (F) HIF1A genes. 

Expression levels were normalized to HPRT1. The differences were evaluated by Mann-Whitney U test. * p < 0.05; ** p < 0.01 and *** p < 0.001 
were considered significant. 



CDSS) and genes involved in the response to hormone 
stimulus (PTGS2, LDLR, PGF, PDGFA, PTGS1, NOS3, 
ADAM9, EGR2, SOCS2, SOCS3, RXRA, ADIPOR1, GAL, 
JUNB, C1QB, RETN, ADM, SORT1, CAV2, CAV1, TNF, 
GRB2, ERBB3, GNG11, BCL2L1, PIK3R1, IL6, IL1RN, 
STATU STAT3, CDKN1A, BMP7, IGFBP2). 

T2D patients presented several genes that were differen- 
tially expressed compared to T1D patients, including genes 
of numerous zinc finger proteins and genes regulating 
transcription (ZNF184, ZNFS76, ZNF449, ZNF594, ZNF641, 
ZNF100, ZNF644, ZNF189, ZNF785, ZNF436, ZNF613, 
among others). Finally, the comparison between T2D and 
GDM showed multiple genes induced in T2D patients, 
including genes involved in DNA-binding, as zinc finger 
(ZNF582, ZNF250, ZNF184, ZNF181, ZNF775, ZNF773, 
ZNF248, ZNF397, ZNF578, ZNF442, ZNF443, among 
others). 

Comprehensive functional analysis using a module 
map approach 

To identify the influence of patient features (array or ex- 
perimental sets) on gene information (gene sets), we 
constructed several module maps, stratifying patients 
according to demographic, clinical, laboratory and 
therapeutic characteristics. We used specialized databases 
associated with diabetes complications [42], gene clusters 
associated with diabetes obtained from association studies 
(GWAS) [43], and isolated immune cell types associated 



with the pathogenesis of diabetes [44]. The most relevant 
modules are shown in Figure 6. 

GDM patients exhibited up-regulated genes observed 
in diabetes complications (including angiopathy) and in 
macrophages. GDM, number of gestations per patient and 
gestation time were associated with the induction of dia- 
betic complications genes. Interestingly, a history of 2 or 
more gestations was positively associated with the modu- 
lation of T2D genes, while in T1D patients the number of 
gestations did no influence the transcription profile of fe- 
male patients (Phenopedia). T1D patients exhibited in- 
duced genes typical of those displayed by dendritic cells 
and repressed genes typical of those presented by B- 
lymphocytes. In addition, T1D patients exhibited induc- 
tion of genes related to diabetes nephropathy. The use of 
insulin did not influence gene expression patterns; how- 
ever, increased serum glucose level was associated with 
the induction of genes related to diabetic retinopathy. In 
patients with T2D, the disease itself as well as the use of 
metformin was associated with the repression of genes im- 
plicated in obesity and diabetic complications. In contrast, 
genes found by T2D GWAS were associated with genes 
induced in T2D patients and the use of metformin. Male 
patients who had T2D for 11 or more years were at 
increased risk of neuropathy. The heatmap shown in 
Figure 6A illustrates these results. 

In the paired analysis of T1D and GDM patients, both 
obesity and a history of two or more gestations were 
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positively associated with T2D development. Conversely, 
in T1D patients, these variables were negatively associ- 
ated. Additionally, in T1D patients, the use of insulin re- 
pressed genes of angiopathy and retinopathy, and 
macrophage gene expression was associated with GDM 
(Figure 6B). In the paired analysis of T2D and GDM, 
the number of gestations was positively associated with 
genes related to diabetic complications. The use of 
metformin was negatively associated with complica- 
tions, and these same genes were repressed in T2D 
(Figure 6C). The paired analysis of T1D and T2D pa- 
tients showed repression of dendritic cell genes associ- 
ated with T2D and the induction of dendritic cell 
genes associated with time of disease. Obesity genes 
were also induced in T1D patients (Figure 6D). 

Confirmation by real-time PCR 

We performed the confirmation of EGR2, RGS1, FOX03A, 
HIF1A, IL1B and SOD2 genes (Figure 7), which were 
chosen because some of these genes have been previously 
described as specific for GDM compared with T2D 
(SOD2, FOX03A, HIF1A). Compared to GDM, these 
genes were downregulated in T2D as evaluated by micro- 
arrays and qRT-PCR, presenting fold changes of 4.34, 1.75 
and 1.72, respectively (microarrays), and 1.51, 2.66 and 
2.31, respectively (qRT-PCR). In addition, compared to 
T2D, the EGR2, IL1B and RGS1 genes were downregu- 
lated when compared to T1D (fold changes of 1.82, 13.96 
and 3.08, respectively, for microarrays and 2.59, 47.83, and 
9.95, respectively, for qRT-PCR), and downregulated when 
compared to GDM (fold changes of 2.22, 11.11 and 4.34, 
respectively, for microarrays, and 2.82, 14.50, and 18.36, 
respectively, for qRT-PCR). 

Discussion 

Diabetes mellitus is one of the most studied diseases, and 
a large amount of information is available in the public da- 
tabases regarding genetic association, meta-analysis and 
associated complications. Few studies, however, have 
systematically compared the major types of diabetes in 
terms of gene expression profiles at the genomic level. 
In this study, we performed an integrative analysis of 
the transcriptome profiles of the major types of diabetes 
using several bioinformatics tools. The major innovation 
of this study was the construction of informative mod- 
ule maps integrating epidemiologic, clinical, laboratory, 
pathogenetic, genetic, and molecular factors implicated 
in diabetes to identify individual and shared features in 
T1D, T2D and GDM. 

The major finding of our global partitioning analysis 
that examined the differences in transcription among each 
diabetes patient was the cluster of genes associated with 
inflammation. The high expression levels of these genes in 
some T1D and GDM patients seem to influence the global 



gene expression pattern of diabetic patients. Indeed, 
several important molecular mechanisms identified by 
clustering account for an intricate array of inflammation 
pathways. One important finding was the up-regulation 
of many mediators of the NOD-like receptor signaling 
pathway. It has recently been suggested that NOD-like 
receptors can be induced by hyperglycemia and oxidative 
stress products, which could link metabolism and inflam- 
mation, particularly through the participation of IL1B 
[45]. In this study, many genes involved in these pathways 
were induced in T1D and GDM patients, including 
NLRP3 (an important receptor of the NOD-like pathway), 
IL1B, CXCL1, CXCL2, IL6, IL8, TNF, RIPK2, TNFAIP3 and 
NFKB1A. TNF and IL1B are strongly involved in the 
regulation of nitric oxide biosynthetic processes [46] 
and are regulated by SOD2 [47]. In addition, other up- 
regulated genes including chemokines (CCL3 } CCL4) } 
cytokines (IL6, IL8, TNF) and transcription factors 
(NFKBIA, MAPK8) are involved with other inflammatory 
processes, such as toll-like receptors signaling pathway, 
and the expression of IL1B and SOD2 genes was con- 
firmed by real-time PCR. Curiously, genes associated with 
these processes were downregulated in T2D compared 
with the other types of diabetes. To understand this result, 
we took advantage of the rank product analysis (paired 
analyses), which showed particular features of the involve- 
ment of these inflammatory pathways in each type of dia- 
betes. The comparisons between T2D and GDM as well 
as between T1D and T2D revealed several interesting re- 
sults. At the same time that NOD-like receptors were in- 
duced in T1D and GDM, hundreds of zinc finger protein 
genes were induced in T2D. Additionally, drug treatment 
with metformin in T2D patients seems to influence the 
gene expression patterns, whereas insulin treatment did 
not. As literature findings have indicated that the expres- 
sion of transcription factors associated with T2D can be 
induced by hypoglycemiants, it is possible to hypothesize 
that these drugs may induce alterations in the normal 
T2D expression profile. Finally, the comparisons between 
T2D and T1D reveal the repression of important genes as- 
sociated with T2D diabetes, particularly transcription fac- 
tors involved in glucose homeostasis (TCF7L2) [24,48], 
NAD + (FOX03) [49] and regulation of cellular and sys- 
temic response to hypoxia (HIF1A) [50], and the differen- 
tial expression of HIF1A and FOX03A in our series was 
also confirmed by real-time PCR. 

To further understand the close similarity between the 
transcription profiles in T1D and GDM, we examined 
the paired comparisons between these types. T1D pa- 
tients exhibited modulation of genes in the MHC region, 
including the induction of HLA-DQA1 and HLA-DQA2. 
The DQA gene encodes the alpha chain of the HLA-DQ 
heterodimeric molecule. HLA-DQ is non-covalently as- 
sociated with DQp, which has been associated with T1D 
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susceptibility in multiple populations [51-54]. In con- 
trast, few studies have associated any HLA class II genes, 
including HLA-DRB3 genes, induced in GDM [53]. Al- 
though a role of HLA class II genes in T1D pathogenesis 
has not yet been established, decreased surface expres- 
sion of HLA-DQ molecules on CD4 and CD8 peripheral 
cells has been reported in recently diagnosed T1D pa- 
tients exhibiting DQB1 susceptibility alleles. This finding 
was attributed to the instability of these molecules on 
the cell surface [54]. As none of our patients were re- 
cently diagnosed, it is interesting to observe that the ex- 
pression of DQA MHC molecules is still modulated 
even after the disease has been present for long periods 
of time. On the contrary, there is also evidence that 
humoral responses to autoantigens may be driven by 
the HLA-DQA1 genes even in recently diagnosed patients 
[55]. MHC susceptibility alleles associated with T1D 
(DQArOS:01, DQB 1*02:0 1/DQB 1*03:02) are different 
from those associated with GDM [56]. In addition to 
MHC genes, T1D patients also showed different levels 
of the killer immunoglobulin-like receptor (KIR3DL2 
and K1R2DS4) family genes compared with GDM patients 
(K1R2DL4). There are no gene expression studies implicat- 
ing KIRs in diabetes. These findings suggest that the close 
similarity of the T1D and GDM transcription profiles may 
be due to the overall inflammatory gene patterns ob- 
served in both conditions, including the modulation of 
several genes primarily involved with the innate im- 
mune response. 

Previously, the construction of module maps for can- 
cer patients has revealed important modules that 
characterize different cancer lineages [41], allowing the 
identification of cancer biomarkers. Major findings re- 
garding each group of diabetes patients are discussed 
below, including epidemiological, clinical, laboratory, 
genetic and pathogenetic features. In all module maps in 
which the gene profile of GDM was compared with 
other types of diabetes, we observed an up-regulation of 
genes typically expressed by macrophages. Overall, the 
comparison with genesets with arraysets yielded positive 
associations, except for T1D, which may be attributed to 
differences in patient age and the autoimmune nature of 
T1D. These macrophage genes were co-regulated with 
those appearing in diabetic complications such as angio- 
pathy and retinopathy, including the 1L1B and RGS1 
genes, as seen in Figure 6A, B and C. As we also observed 
that increased glucose levels were associated with the de- 
velopment of retinopathy, it is possible that the effect of 
high glucose levels on macrophages might be involved in 
the creation of GDM complications. Moreover, GDM 
patients exhibited gene profiles similar to those reported 
for obesity (data available in the Phenopedia public data 
banks), and the majority of patients in our study were 
overweight. It is also interesting to observe that obesity, 



two or more gestations per patient, and gestations 
lasting over 30 weeks exhibited the same modules 
of induced genes in GDM patients, as observed in 
Figure 6A. Additionally, in GDM patients, having two 
or more gestations was positively associated with the 
development of T2D. Indeed, it is important to note 
that obese GDM patients have an increased risk of deve- 
loping T2D [57]. This information is important for future 
studies of GDM. 

In T1D patients, we observed a positive association with 
the profile reported for dendritic cells (ImmGen) 
(Figure 6A), particularly in patients exhibiting long-term 
disease (Figure 6D). As dendritic cells play an important 
role in antigen presentation via MHC class II molecules 
and as MHC class II genes are also induced in long-term 
T1D patients, it is possible to hypothesize that abnormal 
antigen presentation (foreign or self) is a chronic 
phenomenon in T1D. In contrast, T1D patients also ex- 
hibited a repression of genes associated with Bla and Bib 
lymphocytes (naive B cells), in which the Bla lymphocyte 
subtype (CD5+) has been associated with the production of 
natural and autoantibodies [58]. Unfortunately, until the 
moment, there is no public data regarding the gene pattern 
of activated B cells. Another important finding is the simi- 
larity of T1D gene expression profiles with that observed 
for diabetes nephropathy (Phenopedia). As no patients in 
our analysis exhibited clinical nephropathy, it may be valu- 
able to further study the particular patients exhibiting ne- 
phropathy related gene patterns. 

Finally, T2D patients exhibited gene expression pro- 
files that were in disagreement with those reported in 
the public databases (Phenopedia) for diabetic complica- 
tions (Figure 6C). Considering that half of our T2D patients 
exhibited higher median glucose levels than our T1D 
and GDM patients, one could expect to find the induc- 
tion of genes associated with diabetic complications. As 
the T2D patients were treated with many types of medi- 
cations besides hypoglycemiants, it is possible that this 
intriguing and unforeseen finding may be a consequence 
of these treatments. Another possible explanation may 
be related to inflammation genes, which were down- 
regulated in T2D compared to T1D and GDM patients. 
The coupling of inflammation and drugs used to treat T2D 
might be responsible based on these findings. Indeed, the 
use of metformin and other medications (aspirin, captopril, 
atorvastatin, and hydrochlorothiazide) seems to modulate 
the expression of a large number of genes (Figure 6C 
and D), possibly affecting the inflammation status of T2D 
patients. Some of the down- regulated genes observed in 
T2D patients being treated with several drugs included 
1L1B, 1L4, 1L8, CCL2 and INF. All of these genes are in- 
volved in inflammation and are also modulated in macro- 
phages (data not shown). However, we cannot disregard 
the participation of the inflammatory pathway in T2D, as 
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genes that participate in the NOD-like receptors signaling 
pathway probably also play a role in T2D [45]. Recent 
findings have shown that the down- regulation of FOXOl 
expression in macrophages blocks lipid accumulation in 
these cells, affecting many processes [59]. In our analysis, 
FOX03A was also down- regulated in T2D patients com- 
pared to GDM. Both FOX01A and FOX03A genes are 
mediators of the same signaling pathway, AMPK, which 
can be activated as a consequence of long-term metformin 
use [49], and the expression of FOX03A was confirmed 
by real-time PCR. 

In conclusion, our analysis revealed that T1D and 
GDM exhibited a similar up-regulation of inflamma- 
tory genes. 

Methods 

Study population 

We studied 56 adult diabetic patients, 19 presenting T1D 
(7 women/ 12 men) with ages of 18-36 years, 20 presenting 
T2D (13 women/7 men) with ages of 41-72 years and 17 
presenting gestational diabetes with ages of 23 to 40 years. 
For our analysis, we used the mean age of the three groups 
(37 years with SD ± 14). The mean length of disease was 
11 ±5.3 years, while the mean values of glucose and 
HblAc levels were 130.7 mg/dL with SD ± 75.4 and 8.9 
with SD ± 1.8, respectively. For gestational diabetes pa- 
tients, the mean period of gestation was 30.5 with SD ± 
5.7 weeks, and the mean number of pregnancies was 2 ± 1 
per patient. Twelve GDM women had more than two ges- 
tations, while 9 T2D and none T1D women presented this 
characteristic. T1D and GDM patients were treated only 
with insulin, while T2D patients were treated with insulin 
in combination with metformin, captopril, aspirin, atorva- 
statin and hydrochlorothiazide. All patients underwent 
follow-up examinations at the Outpatient Clinics of the 
Division of Endocrinology, Faculty of Medicine of Ribeirao 
Preto, University of Sao Paulo, Brazil. The exclusion 
criteria were based on recent episodes of ketoacidosis, 
active nephropathy, proliferative retinopathy, diabetic 
foot syndrome, high LDL levels and diagnosed cardio- 
vascular diseases. Figure 1 shows a schematic heatmap 
with all demographic, clinical and laboratory patient 
features. The study protocol was approved by the local 
Ethics Committee (Comite de Etica em Pesquisa do Hos- 
pital das Clinicas e da Faculdade de Medicina de Ribeirao 
Preto da Universidade de Sao Paulo, under the permit 
# 9153/2008), and informed consent was obtained from 
all participants. 

Blood collection, peripheral mononuclear cell isolation 
and RNA extraction 

A total of 20 mL of peripheral blood was collected and 
used for the isolation of PBMCs by discontinuous gradi- 
ent density centrifugation on a Ficoll-Hypaque cushion 



(Sigma, St. Louis, MO). Total RNA was extracted using 
the Trizol reagent (Invitrogen, Carlsbad, CA) according 
to the manufacturers instructions. RNA concentrations 
and ratios were checked using a NanoDrop ND-1000 
spectrophotometer (NanoDrop Products, Wilmington, 
DE), and the RNA integrity was assessed by microfluidic 
electrophoresis using a 2100 Bioanalyzer and RNA 6000 
nanochips (Agilent Technologies, Santa Clara, CA). We 
used only samples that exhibited median RNA integrity 
number (RIN) > 9.0. 

Oligo microarrays 

Hybridizations onto whole human genome 4x44K oligo 
microarrays (G4112F, Agilent) were performed using the 
one color (Cy3) Quick Amp labeling kit (Agilent). Briefly, 
500 ng of total PBMC RNA plus spike-in controls were re- 
verse transcribed into double stranded cDNA. The 
primers used for this reaction contained many consecutive 
thymine bases attached to a T7 promoter that paired at 
the 5' end of the first strand of cDNAs. Next, the T7 poly- 
merase was added along with nucleotides labeled with 
fluorescent cyanine-3 (Cy3) dye, which amplified the 
anti-sense complementary RNAs (cRNA). The cRNAs 
were purified and then hybridized to the microarray for 
17 hours at 65°C. After washing, the slides were scanned 
using a DNA Microarray Scanner with Surescan High- 
Resolution Technology (Agilent). A complete file pro- 
viding microarray data from all samples used in this 
study, as well as the numerical quantitative data and 
experimental conditions, is available on line at the 
ArrayExpress public database [60] through the follo- 
wing accession numbers: T1D (E-MEXP-3348), T2D 
(E-MEXP-3287) and GDM (E-MEXP-3349). These data 
correspond to a part of our laboratory databank, which 
were made publicly available and can be used for fur- 
ther studies and/or reanalysis. 

Data quantification and normalization 

Data quantification and quality control were analyzed using 
the Feature Extraction (FE) software version 10.7 (Agi- 
lent Technologies). Expression data were loaded into 
an R-environment [61] using the AgiND package [62], 
a tool developed by the Technologies Avancees pour le 
Genome et la Clinique (INSERM U1090, Marseille, 
France). The AgiND tool is available by request [63]. 
The background adjustment was performed by sub- 
tracting median background values from the median 
expression values obtained by FE and then converting 
the results to log-scale. For each sample, any negative 
values observed were replaced by randomly selected 
small positive values. After log-transformation, quan- 
tile normalization using the normalizeQuantile func- 
tion was performed. 
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Density based filtering and Markov clustering (DBF-MCL) 
and principal component analysis (PCA) 

The Density Based Filtering and Markov Clustering (DBF- 
MCL) algorithm implemented in the Rtools4TB bio- 
conductor package [64,65] was used to extract sets of 
co-regulated genes from our microarray dataset DBF- 
MCL is a tree-step adaptive algorithm that finds gen- 
omic elements and genes located in dense areas, uses 
selected items to construct a graph and finally creates 
a partitioning graph using the Markov Clustering Al- 
gorithm (MCL). A 10% false discovery rate (FDR) was 
used. The principal component analysis (PCA) was 
performed from the informative genes identified by 
DBF-MCL. In this analysis, the R function prcomp [66] 
was used to evaluate PCA and the rgl package [67] to 
construct 3D graphics. 

Statistical analysis 

Differentially expressed genes were identified using the R 
package RankProd [68,69]. Although the non-parametric 
rank product method does not make any assumptions 
about the data distribution, it can provide frequency rank- 
ing scores at each data point and is thus a robust tool for 
creating ranking lists. Genes were considered significantly 
expressed when they presented p-values and percentage of 
false positive predictions (PFP) smaller than 0.001 and 0.05, 
respectively. This test was used to perform paired analysis 
between T1D versus GDM, T2D versus GDM and T1D 
versus T2D. Venn diagrams illustrating up and down- 
regulated genes in each analysis are shown in Figure 4. 

Module-map construction 

Significant and differentially expressed genes obtained 
by the use of DBF-MCL and Rank Products were ini- 
tially clustered using the Cluster 3.0 and TreeView soft- 
wares [70,71]. Then, module maps were constructed 
using an ensemble of tools provided by the Genomica 
software [41], which searches for higher-order modules 
of gene sets and samples. Initially, this algorithm uses 
gene sets to annotate genes that were up-regulated (or 
down-regulated) by at least 1.5-fold. The algorithm com- 
pares the modulated genes with array sets, including 
groups of compartmentalized genes, and organizes them 
into modules (module maps) discriminating variable- 
specific gene patterns according to patient features. It 
was used P value < 0.05 with false discovery ration 
(FDR) < 0.05. 

Module map array set variables 

The variables used to create the experimental sets in- 
cluded demographic parameters (age and gender), clin- 
ical variables (disease duration, obesity, duration and 
number of gestations), laboratory data (serum glucose 
and glycated hemoglobin levels), and treatment features 



(use of hypoglycemiants such as insulin, metformin or 
other drugs). All input variables were transformed into 
binary data (0 or 1), according to the nature of the vari- 
able, i.e., qualitative variables were assigned by the ab- 
sence or presence of the characteristic, and quantitative 
variables were assigned by values below or above the 
mean values. The variables used included age, time of 
disease, gender, serum glucose and glycated hemoglobin 
values (HblAc), use of insulin, use of metformin, use of 
any other type of medication, obesity, gestation time (for 
GDM patients) and number of gestations (for GDM pa- 
tients) (Figure 1). 

Module map gene set variables 

Immune cell specific gene sets were obtained by re- 
analyzing raw data from the ImmGen project [44,72] 
(the detailed procedure is provided as Additional file 1). 
Disease-related gene sets were obtained from GWAS in- 
tegrator [43] and Phenopedia [42]. GWAS integrator is a 
compilation of genes offered by GWAS catalogs and da- 
tabases such as HapMap, SNAP and HugeNavigator 
[73]. All of the information about genes related to T1D 
and T2D were considered for use in this study. Pheno- 
pedia is an available database about genetic association 
studies and meta-analysis summarized in the Human 
Genome Epidemiology (HuGE) encyclopedia [74]. Lists 
of genes can be obtained about specific diseases or re- 
lated complications [75]. The available lists of genes 
related to diabetes mellitus include those associated 
with pre-diabetic status, diabetes, type 1 diabetes, type 
2 diabetes, gestational diabetes, diabetic complications, 
obesity, angiopathy, ketoacidosis, nephropathy, neur- 
opathy and retinopathy. 

Functional analysis 

All clusters of co-regulated genes were analyzed by func- 
tional analysis using the Database for Annotation, 
Visualization, and Integrated Discovery (DAVID) version 
6.7 [76]. This approach was used to identify significant 
biological processes and pathways represented by the 
differentially expressed genes [77]. A biological process 
or pathway was considered to be significant if it con- 
tained a minimum of three genes per category featuring 
score values less than 0.05 including the Benjamini- 
Hochberg correction. Moreover, we selected Kyoto 
Encyclopedia of Genes and Genomes (KEGG) provided 
by DAVID analysis [76] to explore biological pathways. 

Confirmation by real-time PCR 

Relevant and differentially expressed genes obtained 
using the microarray data were confirmed using qRT- 
PCR (Table 1). The sequences of selected genes were re- 
trieved from the NCBI GenBank database [78]. All major 
alternative transcripts were considered to design the 
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Table 1 Primers used for confirmation by qRT-PCR of the I LIB, RGS1, EGR2, FOX03A, SOD2 and HIF1A genes 



Gene 


Accesion numbers 


Primer forward 


Primer reverse 


Size (bp) 


II 1 R 
ILIu 


NM_ 


_000576.2 


D "LLAL/AuMLL I I LLMuuAuAft-j 


c' (ZT(Z ATCC-iTAC kCCTCC ATCC-, 3' 
J -o 1 un 1 1 r\^r\\j\j 1 ULn 1 D 


1 1 1 

I z I 


RGS1 


NM_ 


_002922.3 


5 '-TGGCTGGC^GTG A AG ACTA-3 ' 


5'-GATOTCGAGTGCGGAAGTC-3' 


131 


EGR2 


NM_ 


_000399.3 


5'-GGTGACCATCmCCCAATG-3' 


5 '-TATGGGAGATCCAACGACCT-3' 


120 


F0X03A 


NM_ 


_001 455.3 


5 '-GTGCTA AGCAGGCCTC ATCT-3 ' 


5 '-TTGGC A A AGGG^CTCTG-3 ' 


119 


S0D2 


NM_ 


_00 1024465.1 


5'-GACAAACCTCAGCCCTAACG-3' 


5'-TOGACACCAACAGATGCAG-3' 


124 


HIF1A 


NM_ 


_1 81 054.2 


5'-TCAGCTAmGCGTGTGAGG-3' 


5'-AAAACCATCCAAGGCmCA-3' 


107 


GAPDH 


NM_ 


_002046.3 


5 '-CTCTGCTCCTCCTGTOG AC-3 ' 


5'-ACGACCAAATCCGTOACTC-3' 


112 


HPRT1 


NM_ 


_000 194.2 


5'-GACCAGTCAACAGGGGACAT-3' 


5'-CTGCAI IGI 1 1 IGCCAGTGT-3' 


111 



GAPDH and HPRT1 were selected as endogens. 

primers, and the Primer3 web tool [79] was used to se- 
lect pairs of oligonucleotide primers spanning at the 
exon/exon junction. An optimal melting temperature of 
60°C was standardized for all genes. Transcriptional ex- 
pression levels were determined using SYBR Green (Life 
Technologies), in a StepOne Real Time PCR System 
(Applied Biosystems). The experiments were done in du- 
plicates. The GAPDH and HPRT1 genes were used as 
constitutive genes and the selection of the best genes for 
qRT analysis was done using the Housekeeper software 
[80]. The standard curve for each primer was deter- 
mined and only primers with efficiency upper than 90% 
were used. It was applied the AACT method as de- 
scribed by [81], and the minimum expression value was 
used as calibrator. The graphics were constructed in 
GraphPad Prism version 6.0 and statistical analysis in 
the R statistical environment version 3.0.1. Non- 
parametric Mann- Whitney U test was used, considering 
P value < 0.05 to be significant. 

Conclusion 

The present study revealed that epidemiological, clinical, 
laboratory, immunological, genetic and treatment fea- 
tures influenced the transcriptome profiles in the major 
types of diabetes. In addition, inflammation associated 
with macrophage and dendritic cell function may be 
responsible for clustering GDM and T1D patients to- 
gether, while inflammation in T2D may be influenced by 
drug treatment. 
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